rm(list = ls())

## SET YOUR FILE PATH
setwd("~/Dropbox/Immigration/Analysis/Policy Diffusion/Replication Files/")

library(dplyr)
library(survminer)
library(tidyr)
library(standardize)
library(survival)
library(coxme)
library(frailtySurv)
library(frailtyEM)
library(frailtypack)
library(Hmisc)
library(haven)
library(foreign)
library(coefplot)
library(stargazer)
library(simPH)
library(coefplot)
library(xtable)
library(ggplot2)
library(texreg)
library(leaflet)
library(haven)
library(ggplot2)
library(shinyjs)
library(RColorBrewer)

data <- read_dta("dwrap_index.dta")

data.africa <- subset(data, africa == 1)
data.1980 <- subset(data, year >= 1980)
data.1990 <- subset(data, year >= 1990)
data.indp5 <- subset(data, indp5 == 1)
data.indp10 <- subset(data, indp10 == 1)

set.seed(8675309)

## Set Convergence Criteria
coxph.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000, toler.inf = sqrt(.Machine$double.eps), outer.max = 1000)

coxme.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000)

## TABLE A.36

res1 <- matrix(NA,15,7)

m1b <- coxph(formula = Surv(rest1_start, rest1_end, restriction1) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(reststrata) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1b)
res1[1:2,2] <- c(m1b$coef[1],sqrt(m1b$var[1,1]))
res1[3:4,2] <- c(m1b$coef[2],sqrt(m1b$var[2,2]))
res1[5:6,2] <- c(m1b$coef[3],sqrt(m1b$var[3,3]))
res1[13,2] <- m1b$n
res1[14,2] <- m1b$loglik[2]
res1[15,2] <- AIC(m1b)

m1c <- coxph(formula = Surv(rest1_start, rest1_end, restriction1) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(reststrata) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1c)
res1[1:2,3] <- c(m1c$coef[1],sqrt(m1c$var[1,1]))
res1[3:4,3] <- c(m1c$coef[2],sqrt(m1c$var[2,2]))
res1[5:6,3] <- c(m1c$coef[3],sqrt(m1c$var[3,3]))
res1[13,3] <- m1c$n
res1[14,3] <- m1c$loglik[2]
res1[15,3] <- AIC(m1c)

m1d <- coxph(formula = Surv(rest15_start, rest15_end, restriction_15) ~
               std_ihs_gdppc_lag1 +
                 neggdpshock_lag1+
                 ihs_pop_lag1 +
                 
                 cluster(ccode2) +
                 strata(reststrata15)+
                 frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
               data = data,
               method = c("efron"),
               model = TRUE)

summary(m1d)
res1[1:2,4] <- c(m1d$coef[1],sqrt(m1d$var[1,1]))
res1[3:4,4] <- c(m1d$coef[2],sqrt(m1d$var[2,2]))
res1[5:6,4] <- c(m1d$coef[3],sqrt(m1d$var[3,3]))
res1[13,4] <- m1d$n
res1[14,4] <- m1d$loglik[2]
res1[15,4] <- AIC(m1d)

m1e <- coxph(formula = Surv(rest75_start, rest75_end, restriction_75) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(reststrata75)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1e)
res1[1:2,5] <- c(m1e$coef[1],sqrt(m1e$var[1,1]))
res1[3:4,5] <- c(m1e$coef[2],sqrt(m1e$var[2,2]))
res1[5:6,5] <- c(m1e$coef[3],sqrt(m1e$var[3,3]))
res1[13,5] <- m1e$n
res1[14,5] <- m1e$loglik[2]
res1[15,5] <- AIC(m1e)

m1f <- coxph(formula = Surv(rest5_start, rest5_end, restriction_50) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(reststrata50)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1f)
res1[1:2,6] <- c(m1f$coef[1],sqrt(m1f$var[1,1]))
res1[3:4,6] <- c(m1f$coef[2],sqrt(m1f$var[2,2]))
res1[5:6,6] <- c(m1f$coef[3],sqrt(m1f$var[3,3]))
res1[13,6] <- m1f$n
res1[14,6] <- m1f$loglik[2]
res1[15,6] <- AIC(m1f)

m1g <- coxph(formula = Surv(rest25_start, rest25_end, restriction_25) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(reststrata25)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1g)
res1[1:2,7] <- c(m1g$coef[1],sqrt(m1g$var[1,1]))
res1[3:4,7] <- c(m1g$coef[2],sqrt(m1g$var[2,2]))
res1[5:6,7] <- c(m1g$coef[3],sqrt(m1g$var[3,3]))
res1[13,7] <- m1g$n
res1[14,7] <- m1g$loglik[2]
res1[15,7] <- AIC(m1g)

res1 <- round(res1,3)

latex(res1, file = "")
summary(m1b)
summary(m1c)
summary(m1d)
summary(m1e)
summary(m1f)
summary(m1g)

## TABLE A.37

res2 <- matrix(NA,15,7)

m2b <- coxph(formula = Surv(ewrest1_start, ewrest1_end, ew_restriction1) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(ewreststrata) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2b)
res2[1:2,2] <- c(m2b$coef[1],sqrt(m2b$var[1,1]))
res2[3:4,2] <- c(m2b$coef[2],sqrt(m2b$var[2,2]))
res2[5:6,2] <- c(m2b$coef[3],sqrt(m2b$var[3,3]))
res2[13,2] <- m2b$n
res2[14,2] <- m2b$loglik[2]
res2[15,2] <- AIC(m2b)

m2c <- coxph(formula = Surv(ewrest1_start, ewrest1_end, ew_restriction1) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(ewreststrata) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2c)
res2[1:2,3] <- c(m2c$coef[1],sqrt(m2c$var[1,1]))
res2[3:4,3] <- c(m2c$coef[2],sqrt(m2c$var[2,2]))
res2[5:6,3] <- c(m2c$coef[3],sqrt(m2c$var[3,3]))
res2[13,3] <- m2c$n
res2[14,3] <- m2c$loglik[2]
res2[15,3] <- AIC(m2c)

m2d <- coxph(formula = Surv(ewrest15_start, ewrest15_end, ew_restriction_15) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(ewreststrata15)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2d)
res2[1:2,4] <- c(m2d$coef[1],sqrt(m2d$var[1,1]))
res2[3:4,4] <- c(m2d$coef[2],sqrt(m2d$var[2,2]))
res2[5:6,4] <- c(m2d$coef[3],sqrt(m2d$var[3,3]))
res2[13,4] <- m2d$n
res2[14,4] <- m2d$loglik[2]
res2[15,4] <- AIC(m2d)

m2e <- coxph(formula = Surv(ewrest75_start, ewrest75_end, ew_restriction_75) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(ewreststrata75)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2e)
res2[1:2,5] <- c(m2e$coef[1],sqrt(m2e$var[1,1]))
res2[3:4,5] <- c(m2e$coef[2],sqrt(m2e$var[2,2]))
res2[5:6,5] <- c(m2e$coef[3],sqrt(m2e$var[3,3]))
res2[13,5] <- m2e$n
res2[14,5] <- m2e$loglik[2]
res2[15,5] <- AIC(m2e)

m2f <- coxph(formula = Surv(ewrest5_start, ewrest5_end, ew_restriction_5) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(ewreststrata50)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2f)
res2[1:2,6] <- c(m2f$coef[1],sqrt(m2f$var[1,1]))
res2[3:4,6] <- c(m2f$coef[2],sqrt(m2f$var[2,2]))
res2[5:6,6] <- c(m2f$coef[3],sqrt(m2f$var[3,3]))
res2[13,6] <- m2f$n
res2[14,6] <- m2f$loglik[2]
res2[15,6] <- AIC(m2f)

m2g <- coxph(formula = Surv(ewrest25_start, ewrest25_end, ew_restriction_25) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(ewreststrata25)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2g)
res2[1:2,7] <- c(m2g$coef[1],sqrt(m2g$var[1,1]))
res2[3:4,7] <- c(m2g$coef[2],sqrt(m2g$var[2,2]))
res2[5:6,7] <- c(m2g$coef[3],sqrt(m2g$var[3,3]))
res2[13,7] <- m2g$n
res2[14,7] <- m2g$loglik[2]
res2[15,7] <- AIC(m2g)

res2 <- round(res2,3)

latex(res2, file = "")
summary(m2b)
summary(m2c)
summary(m2d)
summary(m2e)
summary(m2f)
summary(m2g)

## FIGURE A.38

armenia <- subset(data, country!="Armenia")
azerbaijan <- subset(data, country!="Azerbaijan")
kazakhstan <- subset(data, country!="Kazakhstan")
kenya <- subset(data, country!="Kenya")
nigeria <- subset(data, country!="Nigeria")
turkey <- subset(data, country!="Turkey")

armenia <- coxph(formula = Surv(rest1_start, rest1_end, restriction1) ~
               std_ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               
               cluster(ccode2) +
               strata(reststrata) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = armenia,
             method = c("efron"),
             model = TRUE)

summary(armenia)

azerbaijan <- coxph(formula = Surv(rest1_start, rest1_end, restriction1) ~
                   std_ihs_gdppc_lag1 +
                   neggdpshock_lag1+
                   ihs_pop_lag1 +
                   
                   cluster(ccode2) +
                   strata(reststrata) +
                   frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
                   frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
                 data = azerbaijan,
                 method = c("efron"),
                 model = TRUE)

summary(azerbaijan)

kazakhstan <- coxph(formula = Surv(rest1_start, rest1_end, restriction1) ~
                      std_ihs_gdppc_lag1 +
                      neggdpshock_lag1+
                      ihs_pop_lag1 +
                      
                      cluster(ccode2) +
                      strata(reststrata) +
                      frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
                      frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
                    data = kazakhstan,
                    method = c("efron"),
                    model = TRUE)

summary(kazakhstan)

kenya <- coxph(formula = Surv(rest1_start, rest1_end, restriction1) ~
                 std_ihs_gdppc_lag1 +
                   neggdpshock_lag1+
                   ihs_pop_lag1 +
                      
                      cluster(ccode2) +
                      strata(reststrata) +
                      frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
                      frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
                    data = kenya,
                    method = c("efron"),
                    model = TRUE)

summary(kenya)

nigeria <- coxph(formula = Surv(rest1_start, rest1_end, restriction1) ~
                   std_ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                 ihs_pop_lag1 +
                 
                 cluster(ccode2) +
                 strata(reststrata) +
                 frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
                 frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
               data = nigeria,
               method = c("efron"),
               model = TRUE)

summary(nigeria)

turkey <- coxph(formula = Surv(rest1_start, rest1_end, restriction1) ~
                  std_ihs_gdppc_lag1 +
                  neggdpshock_lag1+
                   ihs_pop_lag1 +
                   
                   cluster(ccode2) +
                   strata(reststrata) +
                   frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
                   frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
                 data = turkey,
                 method = c("efron"),
                 model = TRUE)

summary(turkey)

outlierplot2 <- multiplot(armenia, azerbaijan, kazakhstan, kenya, nigeria, turkey, secret.weapon=T, innerCI = 2, coefficients = c("std_ihs_gdppc_lag1"), sort = c("magnitude"), names=c("Armenia", "Azerbaijan", "Kazakhstan", "Kenya", "Nigeria", "Turkey"), horizontal = F, color = c("black", "black" , "black" , "black" , "black"), zeroColor ="red", textAngle = 0, title= "", xlab = "Standardized Coefficient  \n", ylab="")
outlierplot2

